library(rdd)
library(texreg)
library(AER)
library(pBrackets)
library(stargazer)
library(xtable)
library(pryr)

the_prefix <- ""

load(
    file=paste0(
        the_prefix,"data/kff_main_replication_data.RData"
        )
)

sdat$AGEN <- as.numeric(as.character(sdat$AGE))

#### add weights
## code to find weights across surveys: health_survey_load_weights.R
## load("kff_survey_weights_Feb2009_to_Aug2017.RData")

## kff_weights$NUMBER <- kff_weights$survey

## sdat <- merge(
##     sdat,
##     kff_weights,
##     by=c("PSRAID","NUMBER"),
##     all.x=T
## )


lout.pre <- lm(FAVOR ~ I(AGEN<65)+AGEN+I(DATN/100)+I((DATN/100)^2)+EDUC+BLACK+HISP+ASIAN+MALE+I(INCOME/1000)+as.factor(PID5),data=sdat[sdat$DATN < 61 & sdat$AGEN > 61 & sdat$AGEN < 68,], weight=WEIGHT)
lout.post <- lm(FAVOR ~ I(AGEN<65)+AGEN+I(DATN/100)+I((DATN/100)^2)+EDUC+BLACK+HISP+ASIAN+MALE+I(INCOME/1000)+as.factor(PID5),data=sdat[sdat$DATN > 60 & sdat$AGEN > 61 & sdat$AGEN < 68,], weight=WEIGHT)
texreg(
    list(lout.pre,lout.post),digits=2,stars=0.05,custom.model.names=c("Pre-Implementation","Post-Implementation"),
    file=paste0(
        the_prefix, "figs/",
        "tableA15_age_rdd_full_model.tex"
    )
   )

lout.pre_no_controls <- lm(FAVOR ~ I(AGEN<65)+AGEN+I(DATN/100)+I((DATN/100)^2),data=sdat[sdat$DATN < 61 & sdat$AGEN > 61 & sdat$AGEN < 68,], weight=WEIGHT)
lout.post_no_controls <- lm(FAVOR ~ I(AGEN<65)+AGEN+I(DATN/100)+I((DATN/100)^2),data=sdat[sdat$DATN > 60 & sdat$AGEN > 61 & sdat$AGEN < 68,], weight=WEIGHT)

lout.pre_no_party <- lm(FAVOR ~ I(AGEN<65)+AGEN+I(DATN/100)+I((DATN/100)^2)+EDUC+BLACK+HISP+ASIAN+MALE+I(INCOME/1000),data=sdat[sdat$DATN < 61 & sdat$AGEN > 61 & sdat$AGEN < 68,], weight=WEIGHT)
lout.post_no_party <- lm(FAVOR ~ I(AGEN<65)+AGEN+I(DATN/100)+I((DATN/100)^2)+EDUC+BLACK+HISP+ASIAN+MALE+I(INCOME/1000),data=sdat[sdat$DATN > 60 & sdat$AGEN > 61 & sdat$AGEN < 68,], weight=WEIGHT)
texreg(list(lout.pre_no_party,lout.post_no_party),digits=2,stars=0.05,custom.model.names=c("Pre-Implementation","Post-Implementation"))


#### robustness check of discontinuity
vars <- c("EDUC","BLACK","HISP","MALE","PID5","INCOME")
rmat <- matrix(NA,length(vars),3)
for(i in 1:length(vars)){
  txt <- paste("lout.robust <- lm(",vars[i],"~ I(AGEN>64)+AGEN+I(AGEN^2)+I(AGEN^3),data=sdat[sdat$AGEN > 61 & sdat$AGEN < 69,],weight=WEIGHT)",sep="")
  eval(parse(text=txt))
  rmat[i,1] <- summary(lout.robust)$coef[2,1]
  rmat[i,2] <- summary(lout.robust)$coef[2,2]
  rmat[i,3] <- summary(lout.robust)$coef[2,3]
}
rownames(rmat) <- vars
colnames(rmat) <- c("Beta","SE","t-value")
print(xtable(rmat,digits=c(0,3,3,3)),
    file=paste0(
        the_prefix, "figs/",
        "tableA16_rdd_confounding_variables_check.tex"
    )
    )


iout <- IKbandwidth(X=sdat$AGEN,Y=sdat$FAVOR,cutpoint=64.5)
###iout
#[1] 2.417309


sdat$DEM <- with(
    sdat,
    abs(((PID-1)/2)-1)
)

lout.pre.post_party <- lm(DEM ~ I(AGEN<65)*I(DATN>61)+AGEN+I(DATN/100)+I((DATN/100)^2)+EDUC+BLACK+HISP+ASIAN+MALE+I(INCOME/1000),data=sdat[sdat$AGEN > 61 & sdat$AGEN < 68,], weight=WEIGHT)
lout.pre.post_party_control_acafavor <- lm(DEM ~ I(AGEN<65)*I(DATN>61) + FAVOR+AGEN+I(DATN/100)+I((DATN/100)^2)+EDUC+BLACK+HISP+ASIAN+MALE+I(INCOME/1000),data=sdat[sdat$AGEN > 61 & sdat$AGEN < 68,], weight=WEIGHT)

lout.pre.post_ind <- lm(PID==2 ~ I(AGEN<65)*I(DATN>61)+AGEN+I(DATN/100)+I((DATN/100)^2)+EDUC+BLACK+HISP+ASIAN+MALE+I(INCOME/1000),data=sdat[sdat$AGEN > 61 & sdat$AGEN < 68,], weight=WEIGHT)
lout.pre.post_ind_control_acafavor <- lm(PID==2 ~ I(AGEN<65)*I(DATN>61) + FAVOR+AGEN+I(DATN/100)+I((DATN/100)^2)+EDUC+BLACK+HISP+ASIAN+MALE+I(INCOME/1000),data=sdat[sdat$AGEN > 61 & sdat$AGEN < 68,], weight=WEIGHT)

stargazer(
    lout.pre.post_party,
    lout.pre.post_party_control_acafavor,
    lout.pre.post_ind,
    lout.pre.post_ind_control_acafavor,
    omit=c("EDUC","BLACK","HISP","ASIAN","MALE","INCOME"),
    ## type="text",
    omit.stat=c("ser","f","adj.rsq", "rsq"),
    ## table.layout = "ts",
    float=F, digits=2,
    out=paste0(
        the_prefix, "figs/",
        "tableA26_age_rdd_with_party_dv.tex"
    )
)


x1 <- rnorm(10000,mean=summary(lout.pre)$coef[2,1],sd=summary(lout.pre)$coef[2,2])
x2 <- rnorm(10000,mean=summary(lout.post)$coef[2,1],sd=summary(lout.post)$coef[2,2])
thediff <- x2-x1
mean(thediff)

pdf(
    paste0(
        the_prefix,"figs/figure4_rdd-diff-12282018-with-survey-weights.pdf"),
   width=5, height=5
)

par(mar=c(4, 4.2, 4, 3))
plot(x=c(1,2),y=c(summary(lout.pre)$coef[2,1],summary(lout.post)$coef[2,1]),
     ylim=c(-.5,.35),xaxt="n",
	ylab="Effect",pch=16,xlim=c(.5,2.5), main="Under 65 favorability\ncompared to 65 or over",
	xlab="Pre-Implementation       Post-Implementation",cex.lab=1.45, bty="n", yaxt="n", cex=1.5, cex.main=1.5)
lines(x=c(1,1),y=c(summary(lout.pre)$coef[2,1]-1.96*summary(lout.pre)$coef[2,2],summary(lout.pre)$coef[2,1]+1.96*summary(lout.pre)$coef[2,2]), lwd=2)
lines(x=c(2,2),y=c(summary(lout.post)$coef[2,1]-1.96*summary(lout.post)$coef[2,2],summary(lout.post)$coef[2,1]+1.96*summary(lout.post)$coef[2,2]), lwd=2)
brackets(x1=1,x2=2,y1=.15,y2=.15)
num1 <- round(summary(lout.post)$coef[2,1]-summary(lout.pre)$coef[2,1],digits=3)
se1 <- round(sd(thediff),digits=3)
text(x=1.5,y=.29,paste("Diff.=",num1,"\n SE=",se1,sep=""), cex=1.5)

par(new=T)

plot(x=c(1.1,2.1),y=c(summary(lout.pre_no_controls)$coef[2,1],summary(lout.post_no_controls)$coef[2,1]),ylim=c(-.5,.35),xaxt="n",
	ylab="",pch=17,xlim=c(.5,2.5),xlab="",cex.lab=1.5, bty="n", cex=1.3, cex.axis=1.3)

lines(x=c(1.1,1.1),y=c(summary(lout.pre_no_controls)$coef[2,1]-1.96*summary(lout.pre_no_controls)$coef[2,2],summary(lout.pre_no_controls)$coef[2,1]+1.96*summary(lout.pre_no_controls)$coef[2,2]), lwd=2)
lines(x=c(2.1,2.1),y=c(summary(lout.post_no_controls)$coef[2,1]-1.96*summary(lout.post_no_controls)$coef[2,2],summary(lout.post_no_controls)$coef[2,1]+1.96*summary(lout.post_no_controls)$coef[2,2]), lwd=2)

abline(h=0,lty=2, lwd=2)

brackets(x1=1.1,x2=2.1,y1=-.32,y2=-.32,type=1,h=-.05)

num2 <- round(summary(lout.post_no_controls)$coef[2,1]-summary(lout.pre_no_controls)$coef[2,1],digits=3)
thediff2 <- rnorm(10000,mean=summary(lout.pre_no_controls)$coef[2,1],sd=summary(lout.pre_no_controls)$coef[2,2])-rnorm(10000,mean=summary(lout.post_no_controls)$coef[2,1],sd=summary(lout.post_no_controls)$coef[2,2])
se2 <- round(sd(thediff2),digits=3)
text(x=1.6,y=-.44,paste("Diff.=",num2,"\n SE=",se2,sep=""), cex=1.5)

dev.off()







sdat$date <- as.Date(paste(sdat$DATE, "01"), "%b%y %d")

shadowtext <- function(x, y=NULL, labels, col='white', bg='black',
                       theta= seq(0, 2*pi, length.out=50), r=0.2, ... ) {

    xy <- xy.coords(x,y)
    xo <- r*strwidth('A')
    yo <- r*strheight('A')

    # draw background text with small shift in x and y in background colour
    for (i in theta) {
        text( xy$x + cos(i)*xo, xy$y + sin(i)*yo, labels, col=bg, ... )
    }
    # draw actual text in exact xy position in foreground colour
    text(xy$x, xy$y, labels, col=col, ... )
}

plot_and_addLines.pryr %<a-% {
 mod <- lm(
    FAVOR  ~
        I(AGEN < 65)*factor(date),
    data=subset(sdat.rdd, !is.na(AGEN)),
    weight=WEIGHT
)
thedates <- sort(unique(subset(sdat.rdd, !is.na(FAVOR) & !is.na(WEIGHT) & !is.na(AGEN) & !is.na(date))$date))
rdd_resids_64 <- predict(
    mod,
    data.frame(
        AGEN = 64,
        date = thedates
    )
)
rdd_resids_65 <- predict(
    mod,
    data.frame(
        AGEN = 65,
        date = thedates
    )
)

agg_favor <- data.frame(
    date = c(thedates, thedates),
    under65 = c(rep(TRUE, length(thedates)), rep(FALSE, length(thedates))),
    favor_m = c(rdd_resids_64, rdd_resids_65)
)

par(mar=c(5, 5, 4, 1))
## with(
##     subset(agg_favor, under65 & !is.na(under65)),
##     plot(
##         date, favor_m, xlim=c(as.Date(start), as.Date(end)),
##         pch=16, col="purple", bty="n", ylab="Favorability (1 - 4)",
##         xlab="Date -- lines: pre and post implementation averages", main=paste0("ACA favorability", append_title), bty="n",
##         cex.lab=2
##     )
## )
## with(subset(agg_favor, !under65 & !is.na(under65)), points(date, favor_m, col="black", pch=15, bty="n", cex = 0.9))

## abline(v=as.Date("2013-10-01"), lty=3, col="red")
## abline(v=as.Date("2014-03-31"), lty=3, col="red")
## abline(v=as.Date("2014-01-01"), col="red")
## legend("bottomright", legend=c("64", "65"), col=c("purple","black"), lty=1, lwd=c(4, 2), pch=c(16, 15), bty="n")
## mtext("Open Enrollment", side=3, at = as.Date("2014-01-01"), col="red")

mod <- lm(favor_m ~ I(date >= "2014-01-01") * as.integer(under65), data = agg_favor)

pred_df <- data.frame(
        date = rep((as.Date(c(start,"2013-12-01"))), 2),
        under65 = rep(c(TRUE, FALSE), each=2)
    )
pred_df$preds <- predict(
    mod,
    newdata=pred_df
)
##     segments(
##     x0=subset(pred_df, under65)$date[1],
##     x1=subset(pred_df, under65)$date[2],
##     y0=subset(pred_df, under65)$preds[1],
##     y1=subset(pred_df, under65)$preds[2],
##     col="purple", lwd=4
## )
## segments(
##     x0=subset(pred_df, !under65)$date[1],
##     x1=subset(pred_df, !under65)$date[2],
##     y0=subset(pred_df, !under65)$preds[1],
##     y1=subset(pred_df, !under65)$preds[2],
##     col="black", lwd=2
## )


pred_df <- data.frame(
        date = rep((as.Date(c("2014-01-01",end))), 2),
        under65 = rep(c(TRUE, FALSE), each=2)
    )
pred_df$preds <- predict(
    mod,
    newdata=pred_df
)

## segments(
##     x0=subset(pred_df, under65)$date[1],
##     x1=subset(pred_df, under65)$date[2],
##     y0=subset(pred_df, under65)$preds[1],
##     y1=subset(pred_df, under65)$preds[2],
##     col="purple", lwd=4
## )
## segments(
##     x0=subset(pred_df, !under65)$date[1],
##     x1=subset(pred_df, !under65)$date[2],
##     y0=subset(pred_df, !under65)$preds[1],
##     y1=subset(pred_df, !under65)$preds[2],
##     col="black", lwd=2
## )


par(mar=c(5, 5, 5.5, 1))
with(
    subset(agg_favor, under65 & !is.na(under65)),
    plot(
        date, favor_m, xlim=c(as.Date(start), as.Date(end)),
        pch=16, col="purple", bty="n", ylab="Favorability (1 - 4)", xlab="Date -- lines: loess smoother", main=paste0("ACA favorability", append_title), bty="n", type="n",
        cex.lab=2.75, cex.main=2.5, cex.axis=2.75
    )
)


 plx_under <- with(
     subset(agg_favor, under65 & !is.na(under65)),
     predict(
         loess(favor_m ~ as.numeric(date), span=0.75), se=T
     )
 )
  plx_not_under <- with(
     subset(agg_favor, !under65 & !is.na(under65)),
     predict(
         loess(favor_m ~ as.numeric(date), span=0.75), se=T
     )
 )

lines(subset(agg_favor, under65 & !is.na(under65))$date,plx_under$fit, col="purple", lwd=6)
x <- subset(agg_favor, under65 & !is.na(under65))$date
y1 <- plx_under$fit - qt(0.975,plx_under$df)*plx_under$se
y2 <- plx_under$fit + qt(0.975,plx_under$df)*plx_under$se
polygon(c(x, rev(x)), y=c(y1, rev(y2)), lty=2, col=rgb(160, 32, 240, 100, maxColorValue=255), border=NA)


lines(subset(agg_favor, !under65 & !is.na(under65))$date,plx_not_under$fit, col="black", lwd=6)
x <- subset(agg_favor, !under65 & !is.na(under65))$date
y1 <- plx_not_under$fit - qt(0.975,plx_not_under$df)*plx_not_under$se
y2 <- plx_not_under$fit + qt(0.975,plx_not_under$df)*plx_not_under$se
polygon(c(x, rev(x)), y=c(y1, rev(y2)), lty=2, col=rgb(0, 0, 0, 100, maxColorValue=255), border=NA)

abline(v=as.Date("2013-10-01"), lty=3, col="red", lwd=3)
abline(v=as.Date("2014-03-31"), lty=3, col="red", lwd=3)
abline(v=as.Date("2014-01-01"), col="red", lwd=3)
legend("bottomright", legend=c("64", "65"), col=c("purple","black"), lty=1, lwd=c(4), bty="n", cex=3)
mtext("Open Enrollment", side=3, at = as.Date("2014-01-01"), col="red", cex=2)


 par(mar=c(5, 5, 5.5, 1))
with(
    subset(agg_favor, under65 & !is.na(under65)),
    plot(
        date, favor_m, xlim=c(as.Date(start), as.Date(end)),
        pch=16, col="purple", bty="n", ylab="Favorability (1 - 4)", xlab="Date", main=paste0("ACA favorability: 64 (treatment group)", append_title), bty="n", type="n", #\nlines: loess smoother approximating quarterly moving average
        cex.lab=2.5, cex.main=2.25, cex.axis=2.5
    )
)

 plx_under <- with(
     subset(agg_favor, under65 & !is.na(under65)),
     predict(
         loess(favor_m ~ as.numeric(date), span=1/10), se=T
     )
 )
  plx_not_under <- with(
     subset(agg_favor, !under65 & !is.na(under65)),
     predict(
         loess(favor_m ~ as.numeric(date), span=1/10), se=T
     )
 )

lines(subset(agg_favor, under65 & !is.na(under65))$date,plx_under$fit, col="purple", lwd=4)

abline(v=as.Date("2013-10-01"), lty=3, col="red")
abline(v=as.Date("2014-03-31"), lty=3, col="red")
abline(v=as.Date("2014-01-01"), col="red")
legend("bottomright", legend=c("64", "65"), col=c("purple","black"), lty=1, lwd=c(4), bty="n", cex=2)
mtext("Open Enrollment", side=3, at = as.Date("2014-01-01"), col="red", cex=1.8)

lines(subset(agg_favor, under65 & !is.na(under65))$date,plx_under$fit, col="purple", lwd=4)
x <- subset(agg_favor, under65 & !is.na(under65))$date
y1 <- plx_under$fit - qt(0.975,plx_under$df)*plx_under$se
y2 <- plx_under$fit + qt(0.975,plx_under$df)*plx_under$se
polygon(c(x, rev(x)), y=c(y1, rev(y2)), lty=2, col=rgb(160, 32, 240, 100, maxColorValue=255), border=NA)

pred_df <- data.frame(
        date = rep((as.Date(c(start,"2013-12-01"))), 2),
        under65 = rep(c(TRUE, FALSE), each=2)
    )
pred_df$preds <- predict(
    mod,
    newdata=pred_df
)
    segments(
    x0=subset(pred_df, under65)$date[1],
    x1=subset(pred_df, under65)$date[2],
    y0=subset(pred_df, under65)$preds[1],
    y1=subset(pred_df, under65)$preds[2],
    col="purple", lwd=4
)



  shadowtext(
    x=as.Date("2012-12-01"),
    y=subset(pred_df, under65)$preds[1],
    labels = format(round(subset(pred_df, under65)$preds[2], 2), nsmall=2),
    bg = "purple", col = "white", cex=2.5
)


pred_df <- data.frame(
        date = rep((as.Date(c("2014-01-01",end))), 2),
        under65 = rep(c(TRUE, FALSE), each=2)
    )
pred_df$preds <- predict(
    mod,
    newdata=pred_df
)

segments(
    x0=subset(pred_df, under65)$date[1],
    x1=subset(pred_df, under65)$date[2],
    y0=subset(pred_df, under65)$preds[1],
    y1=subset(pred_df, under65)$preds[2],
    col="purple", lwd=4
)

 shadowtext(
    x=as.Date("2014-12-01"),
    y=subset(pred_df, under65)$preds[1],
    labels = format(round(subset(pred_df, under65)$preds[2], 2), nsmall=2),
    bg = "purple", col = "white", cex=2.5
)




  par(mar=c(5, 5, 5.5, 1))
with(
    subset(agg_favor, under65 & !is.na(under65)),
    plot(
        date, favor_m, xlim=c(as.Date(start), as.Date(end)),
        pch=16, col="purple", bty="n", ylab="Favorability (1 - 4)", xlab="Date", main=paste0("ACA favorability: 65 (control group)", append_title), bty="n", type="n", #\nlines: loess smoother approximating quarterly moving average
        cex.lab=2.5, cex.main=2.25, cex.axis=2.5
    )
)

 plx_under <- with(
     subset(agg_favor, under65 & !is.na(under65)),
     predict(
         loess(favor_m ~ as.numeric(date), span=1/10), se=T
     )
 )
  plx_not_under <- with(
     subset(agg_favor, !under65 & !is.na(under65)),
     predict(
         loess(favor_m ~ as.numeric(date), span=1/10), se=T
     )
 )


abline(v=as.Date("2013-10-01"), lty=3, col="red")
abline(v=as.Date("2014-03-31"), lty=3, col="red")
abline(v=as.Date("2014-01-01"), col="red")
mtext("Open Enrollment", side=3, at = as.Date("2014-01-01"), col="red", cex=1.8)

lines(subset(agg_favor, !under65 & !is.na(under65))$date,plx_not_under$fit, col="black", lwd=4)
x <- subset(agg_favor, !under65 & !is.na(under65))$date
y1 <- plx_not_under$fit - qt(0.975,plx_not_under$df)*plx_not_under$se
y2 <- plx_not_under$fit + qt(0.975,plx_not_under$df)*plx_not_under$se
 polygon(c(x, rev(x)), y=c(y1, rev(y2)), lty=2, col=rgb(0, 0, 0, 100, maxColorValue=255), border=NA)

 pred_df <- data.frame(
        date = rep((as.Date(c(start,"2013-12-01"))), 2),
        under65 = rep(c(TRUE, FALSE), each=2)
    )
pred_df$preds <- predict(
    mod,
    newdata=pred_df
)

segments(
    x0=subset(pred_df, !under65)$date[1],
    x1=subset(pred_df, !under65)$date[2],
    y0=subset(pred_df, !under65)$preds[1],
    y1=subset(pred_df, !under65)$preds[2],
    col="black", lwd=2
)

shadowtext(
    x=as.Date("2012-06-01"),
    y=subset(pred_df, !under65)$preds[1],
    labels = format(round(subset(pred_df, !under65)$preds[2], 2), nsmall=2),
    bg = "black", col = "white", cex=2.5
 )







pred_df <- data.frame(
        date = rep((as.Date(c("2014-01-01",end))), 2),
        under65 = rep(c(TRUE, FALSE), each=2)
    )
pred_df$preds <- predict(
    mod,
    newdata=pred_df
)

segments(
    x0=subset(pred_df, !under65)$date[1],
    x1=subset(pred_df, !under65)$date[2],
    y0=subset(pred_df, !under65)$preds[1],
    y1=subset(pred_df, !under65)$preds[2],
    col="black", lwd=2
)

shadowtext(
    x=as.Date("2015-06-01"),
    y=subset(pred_df, !under65)$preds[1],
    labels = format(round(subset(pred_df, !under65)$preds[2], 2), nsmall=2),
    bg = "black", col = "white", cex=2.5
 )




}

## ## to reproduce /exactly/ use WEIGHT_r_and_r
sdat$WEIGHT <- sdat$WEIGHT_r_and_r

pdf(
    paste0(
        the_prefix,"figs/figure3_wh_favor_over_time_rdd_v2_just_64_v_65.pdf"
    ),
    width=18, height=7
)

start <- "2010-01-01"
end <- "2018-01-01"

sdat.rdd <- sdat[sdat$AGEN > 63 & sdat$AGEN < 66 & sdat$date >= start & sdat$date < end,]

append_title <- ""

layout(
    matrix(c(1, 2, 1, 3), 2, 2, byrow=T),
    widths=c(1.75,1), heights=c(1,1)
    )

plot_and_addLines.pryr

dev.off()
